
# Programm for estimations on real data #
#########################################

rm(list=ls(all=TRUE))

# Procedure used to get results as tex table

source("C:\\dossiers\\interactive model\\programs\\Restat_final\\print_results.R")


# Library containing the optimization algorithm

library(limSolve)

# Library on graphics

library(graphics)

# Reading of data

library(foreign)
data <- read.dta("C://dataR-130430.dta")

N <- 148
T <- nrow(data)/148



# Procedure used to find synthetic control weights #
####################################################

wipop <- function(Z0,Z1) {

E <- matrix(1,1,ncol(Z0))
F <- 1

b <- 1

G <- diag(ncol(Z0))
H <- matrix(0,ncol(Z0),1)

res <- lsei(A=Z0, B=Z1, E=E, F=F, G=G, H=H, type=2)
w <- res$X

w
}



# Estimation of parameters using synthetic controls #
#####################################################

coefsSC <- function(y,X,expl) {

yy0  <- matrix(y[(N0*T+1):(N*T),1],T,N-N0)
y0   <- yy0[1:(T0-1),]

yNT  <- matrix(0,N0*T,1)

if (expl==0) {
for (j in 1:N0) {
y1 <- matrix(y[((j-1)*T+1):((j-1)*T+T0-1),1],T0-1,1)
ww <- wipop(y0,y1)
yNT[((j-1)*T+1):(j*T),1] <- yy0%*%ww
}
}

if (expl==1) {

#X0 <- matrix(X[(N0+1):N,],ncol(X),N-N0)
X0 <- t(as.matrix(X[(N0+1):N,]))
Xy0 <- rbind(y0,X0)

for (j in 1:N0) {
y1  <- matrix(y[((j-1)*T+1):((j-1)*T+T0-1),1],T0-1,1)
X1  <- matrix(X[j,],ncol(X),1)
Xy1 <- rbind(y1,X1)

ww <- wipop(Xy0,Xy1)
yNT[((j-1)*T+1):(j*T),1] <- yy0%*%ww
}
}

w0 <- rep(1,N0)%x%c(rep(0,T0-1),rep(1,T-T0+1))

ATTSC <- t(w0)%*%(y[1:(N0*T),1]-yNT)/(N0*(T-T0+1))

ATTSC
}



# Recovery of weights from synthetic controls #
###############################################

weightSC <- function(y) {

yy0  <- matrix(y[(N0*T+1):(N*T),1],T,N-N0)
y0   <- yy0[1:(T0-1),]

yNT  <- matrix(0,N0*T,1)

ww <- matrix(0,N-N0,N0)

for (j in 1:N0) {
y1 <- matrix(y[((j-1)*T+1):((j-1)*T+T0-1),1],T0-1,1)
ww[,j] <- wipop(y0,y1)
}

ww
}



# Application of synthetic controls to permutations #
#####################################################

permutSC <- function(y,S) {

yy <- matrix(y,T,N)
ss <- 1:N

ATTSC <- matrix(0,S,1)

for (s in 1:S) {
sss <- sample(ss,N0,replace=FALSE)

yy0 <- yy[,-sss]
yy1 <- yy[,sss]

y0   <- yy0[1:(T0-1),]
yNT  <- matrix(0,N0*T,1)

for (j in 1:N0) {
y1 <- matrix(yy1[1:(T0-1),j],T0-1,1)
ww <- wipop(y0,y1)
yNT[((j-1)*T+1):(j*T),1] <- yy0%*%ww
}

w0 <- rep(1,N0)%x%c(rep(0,T0-1),rep(1,T-T0+1))

ATTSC[s,1] <- t(w0)%*%(matrix(yy1,N0*T,1)-yNT)/(N0*(T-T0+1))
}

ATTSC
}



# Calcul de l'intervalle de confiance � p% des synthetic controls

CISC <- function(y,X,expl,S) {

cSC <- coefsSC(y,X,expl)

yy  <- matrix(y,T,N)
yya <- yy-as.numeric(cSC)*cbind(rbind(matrix(0,T0-1,N0),matrix(1,T-T0+1,N0)), 
                    matrix(0,T,N-N0))

ss <- 1:N

ATTSC <- matrix(0,S,1)

for (s in 1:S) {
sss <- sample(ss,N0,replace=FALSE)
yy1 <- yya[,sss]
yy0 <- yya[,-sss]
y0  <- yy0[1:(T0-1),]
ww  <- matrix(0,N-N0,N0)
yNT <- matrix(0,N0*T,1)

if (expl==0) { 
for (j in 1:N0) {
y1 <- matrix(yy1[1:(T0-1),j],T0-1,1)
ww[,j] <- as.matrix(wipop(y0,y1),N-N0,1)
yNT[((j-1)*T+1):(j*T),1] <- yy0%*%ww[,j]
 }
}
if (expl==1) {
for (j in 1:N0) {
y1 <- matrix(yy1[1:(T0-1),j],T0-1,1)
X0 <- matrix(X[(N0+1):N,],ncol(X),N-N0)
Xy0 <- rbind(y0,X0)
X1  <- matrix(X[j,],ncol(X),1)
Xy1 <- rbind(y1,X1)
ww[,j] <- as.matrix(wipop(Xy0,Xy1),N-N0,1)
yNT[((j-1)*T+1):(j*T),1] <- yy0%*%ww[,j]
 }
}
w0 <- rep(1,N0)%x%c(rep(0,T0-1),rep(1,T-T0+1))
ATTSC[s,1] <- t(w0)%*%(matrix(yy1,N0*T,1)-yNT)/(N0*(T-T0+1))
}

a=sort(ATTSC,decreasing=FALSE)

coef  <- matrix(0,2,1)
coef[1,1] <- a[floor(25*S/1000)]
coef[2,1] <- a[floor(975*S/1000)+1]
coef
}



# Estimation algorithm steps for least-squares estimations #
############################################################

# Step for parameters of explanatory variables

coefX <- function(yp,Xp,iXpXp,betaL,betaF) {
betaX <- iXpXp %*% crossprod(Xp,yp-matrix(betaF %*% t(betaL),nrow(yp),1))
betaX
}

# Step for constant, individual and time additive effects
# where there is no explanatory variable

coefCit0 <- function(ym,yi,yt) {
coefC <- as.numeric(ym)
coefi <- yi - coefC
coeft <- yt - coefC
coeff <- list(coefC,coefi,coeft)
names(coeff) <- c("coefC","coefi","coeft")
coeff
}

# Step for time interactive effects
# where there is no explanatory variable

coefF0 <- function(yp,nbf) {
ww <- matrix(0,T,T)
for (i in 1:(nrow(yp)/T)) {
wi <- yp[((i-1)*T+1):(i*T)]
ww <- ww+crossprod(t(wi))
}
oww <- order(eigen(ww)$values,decreasing=TRUE)
oww <- oww[1:nbf]
betaF <- matrix(sqrt(T)*eigen(ww)$vectors[,t(oww)],T,nbf)
if (betaF[1]<0) {betaF=-betaF}
betaF
}

# Step for individual interactive effects
# where there is no explanatory variable

coefL0 <- function(yp,betaF) {
betaL <- matrix(0,N,ncol(betaF))
for (i in 1:N) {
betaL[i,] <- crossprod(yp[((i-1)*T+1):(i*T)],betaF)/T
}
betaL
}

# Step for constant, individual and time additive effects

coefCit <- function(ym,Xm,yi,Xi,yt,Xt,betaX) {
coefC <- as.numeric(ym - Xm%*%betaX)
coefi <- yi - Xi%*%betaX - coefC
coeft <- yt - Xt%*%betaX - coefC
coeff <- list(coefC,coefi,coeft)
names(coeff) <- c("coefC","coefi","coeft")
coeff
}

# Step for time interactive effects

coefF <- function(yp,Xp,betaX,nbf) {
N <- nrow(yp)/T
ww <- matrix(0,T,T)
for (i in 1:N) {
wi <- yp[((i-1)*T+1):(i*T)]-Xp[((i-1)*T+1):(i*T),]%*%betaX
ww <- ww+crossprod(t(wi))
}
oww <- order(eigen(ww)$values,decreasing=TRUE)
oww <- oww[1:nbf]
betaF <- matrix(sqrt(T)*eigen(ww)$vectors[,t(oww)],T,nbf)
if (betaF[1]<0) {betaF=-betaF}
betaF
}

# Step for individual interactive effects

coefL <- function(yp,Xp,betaX,betaF) {
N <- nrow(yp)/T
betaL <- matrix(0,N,ncol(betaF))
for (i in 1:N) {
betaL[i,] <- crossprod(yp[((i-1)*T+1):(i*T)]-
             Xp[((i-1)*T+1):(i*T),]%*%betaX,betaF)/T
}
betaL
}



# Estimation of parameters using Bai (2009) #
#############################################

coefsbai <- function(y,X,nbf,initL,initF,iter,expl) {

# Quantities computed once and for all

N <- nrow(initL)
T <- nrow(initF)

numi <- seq(1:N) %x% rep(1,T)
numt <- rep(1,N) %x% seq(1:T)
ym   <- apply(y,2,mean)
yi   <- as.matrix(tapply(y,numi,mean),N,1)
yt   <- as.matrix(tapply(y,numt,mean),T,1)
yp   <- y-as.matrix(yi %x% rep(1,T))-as.matrix(rep(1,N) %x% yt)+ym

if (expl==1) {
Xm   <- t(apply(X,2,mean))
Xi   <- matrix(0,N,ncol(X))
for (k in 1:ncol(X)) { Xi[,k] <- tapply(X[,k],numi,mean) }
Xt   <- matrix(0,T,ncol(X))
for (k in 1:ncol(X)) { Xt[,k] <- tapply(X[,k],numt,mean) }
Xp   <- X-Xi%x%rep(1,T)-rep(1,N)%x%Xt+Xm%x%rep(1,N*T)
iXpXp <- solve(crossprod(Xp))
}

if (expl==0) {
cf     <- coefCit0(ym,yi,yt)
coefF  <- coefF0(yp,nbf)
coefL  <- coefL0(yp,coefF)
}

if (expl==1) {
for (i in 1:iter) {
if (i==1) { coefX <- coefX(yp,Xp,iXpXp,initL,initF) }
else	    { coefX <- coefX(yp,Xp,iXpXp,coefL,coefF) }
cf     <- coefCit(ym,Xm,yi,Xi,yt,Xt,coefX)
coefF  <- coefF(yp,Xp,coefX,nbf)
coefL  <- coefL(yp,Xp,coefX,coefF)
}
}

if (expl==0) {
coeff <- list(cf$coefC,cf$coefi,cf$coeft,coefF,coefL)
names(coeff) <- c("coefC","coefi","coeft","coefF","coefL")
}

if (expl==1) {
coeff <- list(coefX,cf$coefC,cf$coefi,cf$coeft,coefF,coefL)
names(coeff) <- c("coefX","coefC","coefi","coeft","coefF","coefL")
}

coeff
}



# Estimation of standard errors for Bai (2009) #
################################################

stdbai_iid <- function(y,X,coefX,coefC,coefi,coeft,coefL,coefF) {

# Quantities computed once and for all

numi <- seq(1:N) %x% rep(1,T)
numt <- rep(1,N) %x% seq(1:T)
ym   <- apply(y,2,mean)
Xm   <- t(apply(X,2,mean))
ymi  <- as.matrix(tapply(y,numi,mean),N,1)
Xmi  <- matrix(0,N,ncol(X))
for (k in 1:ncol(X)) { Xmi[,k] <- tapply(X[,k],numi,mean) }
ymt  <- as.matrix(tapply(y,numt,mean),T,1)
Xmt  <- matrix(0,T,ncol(X))
for (k in 1:ncol(X)) { Xmt[,k] <- tapply(X[,k],numt,mean) }

yp <- y-as.matrix(ymi %x% rep(1,T))-as.matrix(rep(1,N) %x% ymt)+ym
Xp <- X-Xmi%x%rep(1,T)-rep(1,N)%x%Xmt+Xm%x%rep(1,N*T)

MF <- diag(T)-coefF%*%solve(crossprod(coefF))%*%t(coefF)

llm1 <- solve(crossprod(coefL)/N)

eps  <- matrix(0,N*T,1)
som  <- 0

for (i in 1:N) {
somi <- 0
Xi  <- Xp[((i-1)*T+1):(i*T),]
li  <- matrix(coefL[i,],ncol(coefL),1)

eps[((i-1)*T+1):(i*T),1] <- y[((i-1)*T+1):(i*T),1]-X[((i-1)*T+1):(i*T),]%*%coefX-
				    coefF%*%li-coefC-coefi[i,1]-coeft

for (k in 1:N) {
	Xk  <- Xp[((k-1)*T+1):(k*T),]
	lk <- matrix(coefL[k,],ncol(coefL),1)	
	aik <- as.numeric(t(li)%*%llm1%*%lk)
	somi <- somi + aik*MF%*%Xk/N
}
zi <- MF%*%Xi-somi
som <- som + crossprod(zi)/(N*T)
}
sig2 <- as.numeric(crossprod(eps)/(N*T))

stdbai <- matrix(sqrt(diag(sig2*solve(som))/(N*T)),ncol(som),1)
stdbai
}



# Permutation test for individual effects of the treated to belong #
# to the convex hull of individuals effects of the non treated     #
####################################################################

testpermu <- function(y,X,nbf,iter,expl) {

initL <- matrix(0,N,nbf)

initF <- matrix(0,T,nbf)
for (k in 1:nbf) {
initF[k,k] <- sqrt(T)
}

coefs <- coefsbai(y,X,nbf,initL,initF,iter,expl)

l <- cbind(coefs$coefi,coefs$coefL)

dist <- matrix(0,N0,(N-N0+1))

for (i in 1:N0) {

ll <- rbind(l[i,],l[(N0+1):N,])

for (j in 1:(N-N0+1)) {
l0 <- t(ll[-j,])
l1 <- matrix(t(ll[j,]),nbf+1,1)
ww <- wipop(l0,l1)
dist[i,j] <- sqrt(crossprod(l1-l0%*%ww))
}

}

dist
}



# Estimation of parameters and standard errors using FGLS first diff #
######################################################################

coefs1difFGLS <- function(y,X) {
XX <- cbind(X,rep(1,N)%x%rbind(t(rep(0,T-1)),diag(T-1)))
dX <- matrix(0,N*(T-1),ncol(X)+T-1)
dy <- matrix(0,N*(T-1),1)
for (i in 1:N) {
for (t in 1:(T-1)) {
dy[(i-1)*(T-1)+t,1] <- y[(i-1)*T+t+1]-y[(i-1)*T+t]
dX[(i-1)*(T-1)+t,]  <- XX[(i-1)*T+t+1,]-XX[(i-1)*T+t,]
}
}

idXdX <- solve(crossprod(dX))

D  <- 2*diag(T-1)
for (t in 1:(T-2)) {
D[t,t+1] <- -1
D[t+1,t] <- -1
}

coefX <- solve(crossprod(dX))%*%crossprod(dX,dy)

som <- matrix(0,ncol(X)+T-1,ncol(X)+T-1)
for (i in 1:N) {
som <- som+t(dX[((i-1)*(T-1)+1):(i*(T-1)),])%*%D%*%dX[((i-1)*(T-1)+1):(i*(T-1)),]
}

sig2 	<- as.numeric(crossprod(dy-dX%*%coefX)/(2*N*T-sum(diag(idXdX%*%som))))

# FGLS en supposant que les r�sidus sont iid
# Vm1	<- solve(sig2*D)

# FGLS avec une matrice de r�sidus quelconque

Vm1 <- matrix(0,T-1,T-1)
for (i in 1:N) {
eps <- dy[((i-1)*(T-1)+1):(i*(T-1)),]-dX[((i-1)*(T-1)+1):(i*(T-1)),]%*%coefX
Vm1 <- Vm1+crossprod(t(eps))/N
}
Vm1 <- solve(Vm1)

som1 <- matrix(0,ncol(X)+T-1,ncol(X)+T-1)
som2 <- matrix(0,ncol(X)+T-1,1)

for (i in 1:N) {
som1 <- som1+t(dX[((i-1)*(T-1)+1):(i*(T-1)),])%*%Vm1%*%dX[((i-1)*(T-1)+1):(i*(T-1)),]
som2 <- som2+t(dX[((i-1)*(T-1)+1):(i*(T-1)),])%*%Vm1%*%dy[((i-1)*(T-1)+1):(i*(T-1)),]
}

coef <- solve(som1)%*%som2
std  <- as.matrix(sqrt(diag(solve(som1))),nrow(coef),1)

coef <- as.matrix(coef[1:ncol(X),1],ncol(X),1) 
std  <- as.matrix(std[1:ncol(X),1],ncol(X),1)

coeff <- list(coef,std)
names(coeff) <- c("coef","std")

coeff
}



# Estimation using the three procedures #
#########################################

estimation <- function(y,X,N,T,iter,iterEM,nbf) {

cbai      <- matrix(0,2,ncol(X))
cATTCSC   <- matrix(0,2,ncol(X))
cdiff     <- matrix(0,2,ncol(X))

# Valeur initiale des param�tres � estimer

initL <- matrix(0,N,nbf)

initF <- matrix(0,T,nbf)
for (k in 1:nbf) {
initF[k,k] <- sqrt(T)
}

beta  <- coefsbai(y,X,nbf,initL,initF,iter,1)
sbeta	<- stdbai_iid(y,X,beta$coefX,beta$coefC,beta$coefi,beta$coeft,
                        beta$coefL,beta$coefF)

expl <- 1-(ncol(X)==1)
if (expl==0) {
		cbai[1,1] <- beta$coefX 
		cbai[2,1] <- sbeta[1,1]
		}
if (expl==1) { 
		cbai[1,] <- t(beta$coefX) 
		cbai[2,] <- t(sbeta) 
}

# On ne calcule certains estimateurs que pour nbf=3 vu qu'ils ne 
# d�pendent pas de nbf

if (nbf==3) {

if (expl==0) { 
	cATTCSC[1,1] <- coefsSC(y,0,0)
}
if (expl==1) { 
	cATTCSC[1,1] <- coefsSC(y,X[,2:ncol(var$X)],0) 
}

beta  <- coefs1difFGLS(y,X) 
cdiff[1,] <- t(beta$coef[,1])
cdiff[2,] <- t(beta$std[,1])

}

results <- rbind(cbai,cATTCSC,cdiff)

coeff <- list(results)
names(coeff) <- c("results")

coeff
}





# Parameters

iter   <- 20
iterEM <- 1
nbf    <- 1

T0		<- 8

vecT		<- c(rep(0,T0-1),rep(1,T-T0+1))

vecD		<- tapply(data$z,data$i,mean)
vecD		<- as.matrix(as.numeric((vecD>0)),N,1)

N0		<- sum(vecD)

treated	<- order(rep(vecD,each=T),decreasing=TRUE)



########################################
# ESTIMATION SANS VARIABLE EXPLICATIVE #
########################################

X	    	<- as.matrix(vecD %x% vecT)[treated]
X		<- matrix(X)

# EXIT RATE TO A JOB

y		<- as.matrix(log(data$e))[treated]
y		<- matrix(y)

res <- matrix(0,6,5)

for (nbf in 1:5) {
coef <- estimation(y,X,N,T,iter,iterEM,nbf)
res[,nbf] <- coef$results
if (nbf==3) {permutSC <- coef$pSC}
}

res <- round(res*1000)/1000

nl <- c("Interactive effets,","linear",
        "Synthetic controls","",
        "First difference","")

nc <- c("2","3","4","5","6")

lbl <- "Number of factors"

conf <- matrix(0,6,5)

conf[1,] <- res[1,]-qnorm(.975,mean=0,sd=1)*res[2,]
conf[2,] <- res[1,]+qnorm(.975,mean=0,sd=1)*res[2,]

S    <- 10000
expl <- 0

CI <- CISC(y,X,expl,S)

conf[3,3] <- res[3,3]+CI[1,1]
conf[4,3] <- res[3,3]+CI[2,1]

conf[5,] <- res[5,]-qnorm(.975,mean=0,sd=1)*res[6,]
conf[6,] <- res[5,]+qnorm(.975,mean=0,sd=1)*res[6,]

conf <- round(conf*1000)/1000

textable2(x=res,CI=conf,
          file="C:\\table4a_paper",
          label=lbl,nomligne=nl,nomcol=nc,append=FALSE,decim=3,blanc=1)



# EXIT RATE TO UNKNOWN

y		<- as.matrix(log(data$u))[treated]
y		<- matrix(y)

res <- matrix(0,6,5)

for (nbf in 1:5) {
coef <- estimation(y,X,N,T,iter,iterEM,nbf)
res[,nbf] <- coef$results 
if (nbf==3) {permutSC <- coef$pSC}
}

res <- round(res*1000)/1000

nl <- c("Interactive effets,","linear",
        "Synthetic controls","",
        "First difference","")

nc <- c("2","3","4","5","6")

lbl <- "Number of factors"

conf <- matrix(0,6,5)

conf[1,] <- res[1,]-qnorm(.975,mean=0,sd=1)*res[2,]
conf[2,] <- res[1,]+qnorm(.975,mean=0,sd=1)*res[2,]

S    <- 10000
expl <- 0

CI <- CISC(y,X,expl,S)

conf[3,3] <- res[3,3]+CI[1,1]
conf[4,3] <- res[3,3]+CI[2,1]

conf[5,] <- res[5,]-qnorm(.975,mean=0,sd=1)*res[6,]
conf[6,] <- res[5,]+qnorm(.975,mean=0,sd=1)*res[6,]

conf <- round(conf*1000)/1000

textable2(x=res,CI=conf,
          file="C:\\table4b_paper",
          label=lbl,nomligne=nl,nomcol=nc,append=FALSE,decim=3,blanc=1)



# ENTRY RATE

y		<- as.matrix(log(data$entry))[treated]
y		<- matrix(y)

res <- matrix(0,6,5)

for (nbf in 1:5) {
coef <- estimation(y,X,N,T,iter,iterEM,nbf)
res[,nbf] <- coef$results 
if (nbf==3) {permutSC <- coef$pSC}
}

res <- round(res*1000)/1000

nl <- c("Interactive effets,","linear",
        "Synthetic controls","",
        "First difference","")

nc <- c("2","3","4","5","6")

lbl <- "Number of factors"

conf <- matrix(0,6,5)

conf[1,] <- res[1,]-qnorm(.975,mean=0,sd=1)*res[2,]
conf[2,] <- res[1,]+qnorm(.975,mean=0,sd=1)*res[2,]

S    <- 10000
expl <- 0

CI <- CISC(y,X,expl,S)

conf[3,3] <- res[3,3]+CI[1,1]
conf[4,3] <- res[3,3]+CI[2,1]

conf[5,] <- res[5,]-qnorm(.975,mean=0,sd=1)*res[6,]
conf[6,] <- res[5,]+qnorm(.975,mean=0,sd=1)*res[6,]

conf <- round(conf*1000)/1000

textable2(x=res,CI=conf,
          file="C:\\table4c_paper",
          label=lbl,nomligne=nl,nomcol=nc,append=FALSE,decim=3,blanc=1)








##############################################
# ESTIMATION AVEC LE SCORE CROISE AVEC TREND #
##############################################

initL <- matrix(0,N,nbf)

initF <- matrix(0,T,nbf)
for (k in 1:nbf) {
initF[k,k] <- sqrt(T)
}

X	    	<- as.matrix(vecD%x%vecT)[treated]
psorted	<- as.matrix(data$p)[treated]

numi <- rep((1:N),each=T)
psorted	<- matrix(unlist(tapply(psorted,numi,mean)))
XX  <- cbind(X,psorted%x%(1:T)/T)


# EXIT RATE TO A JOB

y		<- as.matrix(log(data$e))[treated]
y		<- matrix(y)

res <- matrix(0,6,5)

for (nbf in 1:5) {

beta  <- coefsbai(y,XX,nbf,initL,initF,iter,1)
sbeta	<- stdbai_iid(y,XX,beta$coefX,beta$coefC,beta$coefi,beta$coeft,
                         beta$coefL,beta$coefF)

res[1,nbf] <- beta$coefX[1,1] 
res[2,nbf] <- sbeta[1,1] 

# On ne calcule certains estimateurs que pour nbf=3 vu qu'ils ne 
# d�pendent pas de nbf

if (nbf==3) {

res[3,nbf] <- coefsSC(y,psorted,1) 

beta  <- coefs1difFGLS(y,XX) 
res[5,nbf] <- beta$coef[1,1]
res[6,nbf] <- beta$std[1,1]
}
}

res <- round(res*1000)/1000

nl <- c("Interactive effets,","linear",
        "Synthetic controls","",
        "First difference","")

nc <- c("2","3","4","5","6")

lbl <- "Number of factors"

conf <- matrix(0,6,5)

conf[1,] <- res[1,]-qnorm(.975,mean=0,sd=1)*res[2,]
conf[2,] <- res[1,]+qnorm(.975,mean=0,sd=1)*res[2,]

S    <- 10000
expl <- 1

CI <- CISC(y,psorted,expl,S)

conf[3,3] <- res[3,3]+CI[1,1]
conf[4,3] <- res[3,3]+CI[2,1]

conf[5,] <- res[5,]-qnorm(.975,mean=0,sd=1)*res[6,]
conf[6,] <- res[5,]+qnorm(.975,mean=0,sd=1)*res[6,]

conf <- round(conf*1000)/1000

textable2(x=res,CI=conf,
          file="C:\\tableC4a_paper",
          label=lbl,nomligne=nl,nomcol=nc,append=FALSE,decim=3,blanc=1)



# EXIT RATE TO UNKNOWN

y		<- as.matrix(log(data$u))[treated]
y		<- matrix(y)

res <- matrix(0,6,5)

for (nbf in 1:5) {

beta  <- coefsbai(y,XX,nbf,initL,initF,iter,1)
sbeta	<- stdbai_iid(y,XX,beta$coefX,beta$coefC,beta$coefi,beta$coeft,
                         beta$coefL,beta$coefF)

res[1,nbf] <- beta$coefX[1,1] 
res[2,nbf] <- sbeta[1,1] 

# On ne calcule certains estimateurs que pour nbf=3 vu qu'ils ne 
# d�pendent pas de nbf

if (nbf==3) {

res[3,nbf] <- coefsSC(y,psorted,1) 

beta  <- coefs1difFGLS(y,XX) 
res[5,nbf] <- beta$coef[1,1]
res[6,nbf] <- beta$std[1,1]
}
}

res <- round(res*1000)/1000

nl <- c("Interactive effets,","linear",
        "Synthetic controls","",
        "First difference","")

nc <- c("2","3","4","5","6")

lbl <- "Number of factors"

conf <- matrix(0,6,5)

conf[1,] <- res[1,]-qnorm(.975,mean=0,sd=1)*res[2,]
conf[2,] <- res[1,]+qnorm(.975,mean=0,sd=1)*res[2,]

S    <- 10000
expl <- 1

CI <- CISC(y,psorted,expl,S)

conf[3,3] <- res[3,3]+CI[1,1]
conf[4,3] <- res[3,3]+CI[2,1]

conf[5,] <- res[5,]-qnorm(.975,mean=0,sd=1)*res[6,]
conf[6,] <- res[5,]+qnorm(.975,mean=0,sd=1)*res[6,]

conf <- round(conf*1000)/1000

textable2(x=res,CI=conf,
          file="C:\\tableC4b_paper",
          label=lbl,nomligne=nl,nomcol=nc,append=FALSE,decim=3,blanc=1)



# ENTRY RATE

y		<- as.matrix(log(data$entry))[treated]
y		<- matrix(y)

res <- matrix(0,6,5)

for (nbf in 1:5) {

beta  <- coefsbai(y,XX,nbf,initL,initF,iter,1)
sbeta	<- stdbai_iid(y,XX,beta$coefX,beta$coefC,beta$coefi,beta$coeft,
                         beta$coefL,beta$coefF)

res[1,nbf] <- beta$coefX[1,1] 
res[2,nbf] <- sbeta[1,1] 

# On ne calcule certains estimateurs que pour nbf=3 vu qu'ils ne 
# d�pendent pas de nbf

if (nbf==3) {

res[3,nbf] <- coefsSC(y,psorted,1) 

beta  <- coefs1difFGLS(y,XX) 
res[5,nbf] <- beta$coef[1,1]
res[6,nbf] <- beta$std[1,1]
}
}

res <- round(res*1000)/1000

nl <- c("Interactive effets,","linear",
        "Synthetic controls","",
        "First difference","")

nc <- c("2","3","4","5","6")

lbl <- "Number of factors"

conf <- matrix(0,6,5)

conf[1,] <- res[1,]-qnorm(.975,mean=0,sd=1)*res[2,]
conf[2,] <- res[1,]+qnorm(.975,mean=0,sd=1)*res[2,]

S    <- 10000
expl <- 1

CI <- CISC(y,psorted,expl,S)

conf[3,3] <- res[3,3]+CI[1,1]
conf[4,3] <- res[3,3]+CI[2,1]

conf[5,] <- res[5,]-qnorm(.975,mean=0,sd=1)*res[6,]
conf[6,] <- res[5,]+qnorm(.975,mean=0,sd=1)*res[6,]

conf <- round(conf*1000)/1000

textable2(x=res,CI=conf,
          file="C:\\tableC4c_paper",
          label=lbl,nomligne=nl,nomcol=nc,append=FALSE,decim=3,blanc=1)











































# GRAPHIC WHEN NBF=1 #
######################

# Horizontal axis: individuel additive effect
# Vertical axis: Multiplicative individual effect

X	    	<- as.matrix(vecD %x% vecT)[treated]
X		<- matrix(X)

# Exit rate to a job

y		<- as.matrix(log(data$e))[treated]
y		<- matrix(y)

initL <- matrix(0,N,nbf)

initF <- matrix(0,T,nbf)
for (k in 1:nbf) {
initF[k,k] <- sqrt(T)
}

expl <- 0
nbf <- 1

coefs <- coefsbai(y,X,nbf,initL,initF,iter,expl)

l <- cbind(coefs$coefi,coefs$coefL)

colorl <- c(rep("red",N0),rep("blue",N-N0))
pchv   <- c(rep(17,N0),rep(0,N-N0))

plot(l[(N0+1):N,1],l[(N0+1):N,2], col="blue", pch=1,
     xlab="Additive local effect", ylab="Multiplicative local effect")

points(l[1:N0,1],l[1:N0,2], col="red", pch=17)





# Permutation test for the individual effects of treated #
# to be in the convex hull of those of non treated	   #
##########################################################

X	    	<- as.matrix(vecD %x% vecT)[treated]
X		<- matrix(X)

# Exit rate to a job

y		<- as.matrix(log(data$e))[treated]
y		<- matrix(y)

nbf <- 2
dist <- testpermu(y,X,nbf,iter,1)

d <- matrix(0,N0,1)

for (i in 1:N0) {
d[i,1]=mean((dist[i,1]>dist[i,2:(N-N0+1)]))
}

ww <- weightSC(y)

wpositif <- matrix(colSums((ww>0)),N0,1)

wd <- cbind(wpositif,d)

wd
